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ABSTRACT 

t Context. To progress in the understanding of the solar wind and coronal dynamics, numerical modeling must include the transition 

region within the simulation domain, and not just as a bottom boundary. Published simulations including a transition region within the 
domain often do not take into account any modeling of the heat sources and sinks which generate the corona; in the rare self-consistent 

CNj simulations including it, the transition region itself appears to be chaotic in space and time. 

Aims. Our aim is to investigate the response to perturbations of the solar atmosphere including transition region and wind, and more 
O specifically how wave propagation is modified by the presence of heat sources and sinks, in the simple ID, hydrodynamical case, 

including chromosphere and solar wind. 

Methods. We integrate the time-dependent hydrodynamic equations of the solar wind with spherical symmetry, including conduction, 
CO radiative cooling and a prescribed mechanical heat flux. Once a quasi-stationary wind is established, we study the response of the 

system to pressure oscillations at the photospheric boundary. We use transparent boundary conditions 

Results. We find that wavepackets with high enough amplitude propagating upward from the photosphere implode just below the 
transition region. This implosion is due to the radiative cooling term generating pressure holes close to the wave crests of the wave, 
which make the wave collapse. In the case where heat sources and sinks are not present in the equations, the wave remains stable 
^ whatever the initial wave amplitude, which is compatible with published work. 

^ Conclusions. The instabiUty found here is not an instability of the TR itself: on the contrary, instability ceases when the wavepacket 

enters the TR where conduction is able to balance the cooling. However, the TR as a whole can be destabilized by such implosions, 
^ which should be observable when and where the TR is high enough above the optically thick regions. 

Key words, hydrodynamics - instabilities - Methods: numerical - Sun: transition region 
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1. Introduction The transition region itself might be the source of inter- 

esting dynamics still to be understood. Indeed, simulations by 

O. , , . , 1 1 J- 1 • • Suzuki and Inutsuka (2005) (a global 1.5D time-dependent so- 

Among the many choices to be made before beginmng a numer- , j^j j .juci-^ufn ^ 

• . , , -,, ■ 1-1 lar wind study, denoted by SI in the tollowine) and Gudiksen 

Oical study ot the corona-solar wind dynamics, one must decide . /r,r,r>c\ /ot^ i f i » u \ j 

, ^ ■' . , ^ , J ■ , , , J • and Nordlund (2005) (3D volume oi solar atmosphere), consider 

— H where to place the bottom boundary; in particular, should we in- c ^^ ^i. ■ c ^c • » » ^ i. ^- tt. 

17^ , J , J ,j, , , I n . 1 J- sucessrully the issue or seU-consistent coronal heatmg. They 

f"-. elude the dense cold layers below the corona, or not .' Includins u ^t. £ j ^t. * ^t. * • • u ^t. • j 

^, .„ . , . . both find that the transition region is erratic both in time and 

them will increase the numerical cost (as the scale height or the 

cold layers is small), and, as well, resolving the chromospheric ^P^''^- 

. ^ transition region itself will require to increase the spatial and ^h^ second question is to specify the energy equation. The 

><; temporal resolutions there. So, we are tempted not to include the mmimal energy equation mcludes the adiabatic term with y = 

5-H cold layers. The work by Lionello et al. (2001) and by Endeve et 5/3. If we are not interested in the dynamics of the transition re- 

^ al. (2003) shows that interesting results can be obtained in this gion, that is all for the energy equation, and we simply start with 

way, albeit with numerous difficulties. ^ gi'^en atmosphere temperature profile, possibly adding extra 

coupling terms in the energy equation to limit the numerical dif- 

However, not including the cold layers means not includ- fusion of the temperature profile. This is done in De Pontieu et 

ing completely the fi-ansition region, so that we are left with al. (2005) in their work on spicules, and more or less so also 

a coronal bottom boundary, where we have to chose some (arbi- (Velh, personal communication) in Del Zanna et al. (2005) in 

trary) conditions. If we are not especially interested by the time- their work on coronal seismology. 

dependent atmosphere, then this is workable (see Lionello et al. Alternatively, if one includes heat sources and sinks, it is rea- 

2001). But if we are, then the response of the corona to pertur- sonable to expect that the dynamics should be modified, even for 

bations might depend much on these arbitrary boundary condi- small amplitude perturbations. There is one single simulation of 

tions. This might show up when perturbing the bottom bound- the corona and solar wind which includes heat sources and sinks 

ary (compare the results obtained by using "line-tied", i.e., fuU self-consistently, the one by SI: it used 14000 grid points, which 

reflexion in Aulanier et al. (2005) and those obtained by using is clearly not generalizable to 2D or 3D. Our aim here is double: 

complete transparency in Grappin et al. (2005)). to devise a reduced version of such a ID model, with a reduced 
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number of grid points, e.g. 300 points instead of 14000, so al- 
lowing future generalization to 2D/3D, and use it to study the 
specific effects on solar atmosphere dynamics (including transi- 
tion region) of heat sources and sinks. 

In this preliminary work, we solve the hydrodynamics equa- 
tions, not the MHD equations, and inject pressure waves at the 
photospheric level. Since we do not expect that plain pressure 
waves are able to heat the corona, we use a prescribed mechani- 
cal energy flux which is a function of heliocentric distance as e.g. 
in Wang, 1994: hence our model is not self-consistent, contrary 
to that of SI. 

A basic starting point for the understanding of waves in the 
solar atmosphere is the work by Fleck and Schmitz (1991) and 
Kalkofen et al. (1994). They studied the response to base pertur- 
bations of an isothermal atmosphere, with only adiabatic terms 
in the energy equation, and a uniform gravity. Recall that ID 
pressure oscillations propagate or not, depending on the fre- 
quency being larger or smaller thant the cut-off Lamb frequency: 



COL = cl(2H) = yglilc) 



(1) 



c being the sound speed, g the gravity, H the pressure scale 
height. When propagative, waves have their energy flux about 
constant in the limit of high frequencies. When evanescent, all 
the atmosphere oscillates in phase. This is the asymptotic sta- 
tionary state. But before this, interesting transients occur, as has 
been discovered by the above authors. Long-Uving transients 
with the cut-off frequency dominate the spectrum in a large 
part of the atmosphere and during a long time, except in the 
very lower layers. Fleck and Schmitz have shown that the phe- 
nomenon is about the same, whether a given transition region 
is present or not. This was proposed to be a possible explana- 
tion for the observed prevalence of 3-min oscillations in chro- 
mospheric lines. How are these results modified when energy 
sources and sinks are taken into account? 

We continue here this work by adding heat sources and sinks 
to the energy equation. We wiU show that the transients with cut- 
off frequency become unstable when they arrive below the TR, 
at an altitude where radiative cooling begins to be active. They 
may implode, leading to catastrophic events. This work will at 
the same time demonstrate the feasibility of a time-dependent 
model of a corona-solar wind starting at the photospheric level 
with a reduced number of grid points. 

The plan is the following. Section 2 describes the equations 
and the method. Section 3 reports the results. The last section is 
a discussion. 



2. Equations and method 

We consider a spherically symmetric solar corona, and solve the 
time-dependent gas equations, which are given by the equation 
for the radial velocity: 

duldt + udu/dr + il/p)dP/dr + GM/r^ = (2) 
Pressure and temperature equations: 

dP/dt + udP/dr + yPdivu = Q (3) 

dTjdt + udT/dr + (y - V)Tdivu = Qjp (4) 

where p - nipn. The density n is deduced from the pressure and 
temperature, using the equation of state: 



T being the average temperature of ions and electrons. The extra, 
non-adiabatic source term Q is given by either one of the two 
following models. 

(a)the "fuU model" including heat sources and sinks: 



e = -(r - ^)idivF„ +p^KT) + divF,) + KiT" 



(6) 



where F^, A(T), and Fc are the (prescribed) mechanical en- 
ergy flux, the radiative cooUng term and the conductive flux, 
and K\T" provides a filtering of the temperature, with T" be- 
ing a second order derivative weighted by the local mesh Ar. 
T" = d^T /dP'iAr/ < Ar >)^, < Ar > being the average mesh, 
(b) the "relaxation model" including as the only non adiabatic 
term a small smoothing term: 



Q = >ciiT"-To{rr) 



(7) 



where T()(r) is a prescribed temperature profile, defining the tran- 
sition region. In both the full and the relaxation models, the co- 
efficient Ki is small enough, so that its characteristic time scale 
is much longer than the time scale of the phenomenon studied. 

The prescribed phenomenological mechanical flux which 
leads to global coronal heating is 



F„(r) = erF°JRJrfexp((R - R,)/Rh) 



(8) 



A typical value for the energy flux is F^ = IQ^erg/cm^/s, and 
the characteristic scale for energy dissipation is Rh = Rg. The 
radiation loss function A(7') is a rough fit to the one by Athay 
(1986). More precisely: 



A = A10-('''«"'<^/^"'V(7') 



(9) 



with Tm = 0.2MK gives the maximum cooling temperature. The 

factor f(T) is between and 1 . It is set to zero if the temperature 
is smaller than a threshold or larger than IMK. Between T^. 
and r** = Q.Q2MK, f(T) varies linearly between and 1. The 
low cut-off temperature is set to: 



= OMMK 



(10) 



Finally, the factor A is fixed so as to lead to a peak value A^ax = 
lO^^^cgs (but see Table 1). The conductive Spitzer-Harm flux is 



F, = -e,KQT^'^dT/dr 



(11) 



where kq = IQ'^cgs. This value is in between the conductivity 
for the proton and electron. This form being highly demanding 
in computer resources (and mesh size), we reduce steepness of 
the resulting flux in two ways. First, we use a linear temperature 
dependence for conductivity at temperature lower than 0.25 MK 
and pass progressively to the Spitzer-Harm 5/2 power law be- 
tween 0.25 MK and IMK (see Linker et al. (2001)). Second, we 
limit the resulting flux by requiring that the associated charac- 
teristic time be not smaller than a prescribed value t/™. The 
final form of the heat flux is thus: 



-1 



(12) 



P = 2nkT 



(5) 



with = K^T'^^^^'^in I Arf- 1 p and Ar is the local mesh size. 
Note that in practice, the limiting factor will only be active in 
the preliminary phase of the runs (run SI in Table 1). 

We use a non-uniform mesh with 300 points following a 
logarithmic progression of ratio q= 1.025. This corresponds to 
a minimum mesh size being Ar = 10"^/?^ = 70 km at the sur- 
face, and a maximum mesh Ar = 0.4/?, at the outer boundary 
(r = l5Rs). Our temporal scheme is Runge-Kutta of order 3; 
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Run 


time 


fO 




KO 


T lim 


N 


SI 


0,20 


4W* 


10-'-' 


10-' 


10-" 


300 


S2 


20,40 


910* 


10-22 


10-*^ 




300 


S3 


40,60 


910" 


10-22 


10-' 


10-^ 


300 


Dl 


60,70 


910-* 


10-22 


10-" 


10-^ 


599 



Table 1. List of runs leading to stationary corona and wind, 
with heat and source coefficients. Runs SI, S2, S3, Dl follow 
each other. SI and S2 respectively correspond to the two phases 
mentioned in the text. S3 is almost identical to S2 (only the vis- 
cous terms differ). Dl starts from an extrapolation of run S3 to 
a more refined mesh and then relaxes. Time: beginning and end 
time of each phase; F° : base mechanical flux; A,„ax' maximum 
of cooling term A(r); /to: conductivity parameter; t/,-,,,: charac- 
teristic time defining the attenuation of the conductive flux; N is 
the resolution 



the spatial scheme is a compact finite difference scheme of or- 
der 6 (Lele, 1992) modified to be able to cope with non-uniform 
grids (J.-M. Le Saout 2003, Grappin et al. 2005). Note however 
that to compute the temperature gradients which appear in the 
conductive term, we use a scheme or order two, which is more 
stable. Finally, we use 6-th order filtering, also defined in Lele 
(1992). Filtering is used in several places. It is used for the veloc- 
ity field (roughly once every 10 time steps), and for smoothing 
different quantities before computing the right-hand sides of the 
equations. Specifically, the (logarithmic) gradients of pressure, 
density and temperature are systematically filtered; finally the 
rhs of the velocity equation is filtered: the latter considerably in- 
creases the stability of the schema. Time-step is automatically 
adapted to the different terms, using dimensional estimations of 
the characteritstic times. 

The boundary conditions are imposed via the characteristic 
form of the equations. They are: no incoming perturbation at the 
inner boundary, and, starting with time t=l, a constant depres- 
sion at the outer boundary, introduced via the ingoing character- 
istics. This depression stops as soon as the sonic Mach number 
is reached, since thereafter no incoming signal can progress into 
the domain from the exterior 

Units used in the following are the solar radius for distance, 
MK for temperature, and km/s or m/s for velocity. Time indi- 
cated in figures and tables is mostly measured in numerical unit, 
which is lh30. The list of runs building the stationary corona and 
wind is given in Table 1, a liste of runs in which we perturb the 
atmosphere with waves injected at the base is given in Table 2. 

3. Results 

3.1. Wave transfer in a relaxation model 

We first consider in this subsection a two-temperature atmo- 
sphere without energy sinks and sources, to check that we re- 
cover basic wave properties in a simple configuration. The atmo- 
sphere is made of two parts: a low, dense part at T = O.OIMK 
and a corona at T - IMK. The transition region lies at about 
R = 0.015/?j. Apart from that, the parameters are the same as 
for the full wind model: domain size L - 14/?s, base density is 
lO'^COT-^*, and same grid. Figure [l] shows the density and tem- 
perature profiles (left). In the right panel, we show the profiles of 
the (square root of the) energy flux, which result from injecting a 
monochromatic wave from time t=0 at the base. The frequency 
is dLt = 100, which is higher than the maximum of cut-off fre- 
quency (about a>o - 90) in the low-temperature region. Hence 
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Runs 


energy 


r 


To 


ev/prop 


amplitude 


IN 


Al 


no (Ij 


6 min 


7.6 min 


prop 


U.4 m/s 


^nn 
jUU 


A 1 ^> 

Ala 


no (ij 


6 min 


7.6 min 


prop 


l.J m/s 


ir\r\ 


A9 


no (,lj 


1 1 min 


7.6 min 


GV 


A m /c 
H- 111/ s 




B4 


yes 


10 min 


6-8 min 


ev 


4m/s 


300 


B40 


yes 


10 min 


6-8 min 


ev 


40m/s 


300 


C40 


no(2) 


10 min 


6-8 min 


ev 


40m/s 


300 


C400 


no(2) 


10 min 


6-8 min 


ev 


400m/s 


300 


D40 


yes 


10 min 


6-8 min 


ev 


40m/s 


599 


D80 


yes 


10 min 


6-8 min 


ev 


80m/s 


599 



Table 2. List of runs with wave injection, runs Al, A2 start in- 
jecting waves at time t=0; other runs start from a relaxed wind 
flow. Runs Bxx, Cxx start from run S3, runs Dxx start from 
run Dl (see above Table 1). Energy: yes means energy equa- 
tion is present; no(l) means coupling of temperature with step- 
wise constant profile; no(2) means coupling temperature with 
the stationary profile resulting from integrating the full energy 
equation, r and tq are resp. the period of the injected wave and 
the cut-off period of the dense layers. Ev/prop: ev means evanes- 
cent in dense layers, prop means propagative in dense layers. 
Amplitude is the base amplitude of the wave. N is the number of 
grid points. 



the quantity plotted should be approximately invariant, except 
for reflexions, due possibly to the frequency being too low in the 
cold region, and most certainly to reflexions at the TR. The the 
sudden energy flux decrease after R - I. IR^ is due to the filter- 
ing, which can be made lighter, so that the damping distance is 
increased. Filtering should certainly be decreased in view of the 
future sef-consistent heating of the corona by wave dissipation. 

run Al run Al 



10"^ 
10" 







1.00 


4x10 ■ 








2x10"' 






0.10 ^ 











-2x10" = 


J 




0.01 


-4x10"' 



z/R_s R-1 

Fig. 1. Waves in a two-temperature atmosphere (relaxation 
model): Left: Density (thick line) and Temperature (plain line) 
profiles; Right: run Al, waves with frequency above cut-off, pro- 
files of (pr^Vg)'^^M in a the (1,1.5) distance range. 



The main property which interests us here is the very dif- 
ferent response of the atmosphere to propagative and evanescent 
waves explained by Fleck and Schmitz (1991) and KaUcofen et 
al. (1994), as mentioned in Section 2. This is shown in Figure|2] 
We compare there the velocity response of the TR to waves re- 
spectively above and below the cut-off (runs Ala and A2). The 
signal at the base is shown as dotted points. One sees that, while 
the oscillation at the TR has the same frequency as injected in the 
high-frequency case, this is not the case in the evanescent case: 
at the TR (and in lower layers as well), the frequency becomes 
very close to the cut-off frequency (period about 7.6 min.). Note 
also that the TR response is largest during the first two or three 
periods. 
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run Ala 



run A2 



40-0 40.2 40-4 40.6 40-8 41. i 
time 




40.0 40-2 40.4 40.6 40.8 41.0 
time 




0.0001 0.0010 0.0100 0.1000 1.0000 10.0000 



Fig. 2. Oscillations at TR with two-temperature atmosphere: ve- 
locity at TR versus time (dotted: base velocity); left: run Ala, 
propagative incident wave with period 6 min. (below cut-off pe- 
riod - 7.6 min.); right: run A2, evanescent waves with period 1 1 
min. 



3.2. Generation of full wind model 

We give here an account of the method used to generate the 
wind and hot corona. We start with an atmosphere uniformly 
at 6000 /T, and a corresponding hydrostatic equilibrium. The ini- 
tial hydrostatic stratification is too strong to start immediately to 
integrate the whole set of fluid equations, essentially because the 
density is very low at the distance of R = 15-Rj, which leads to 
very small characteristic times. We thus decide to first raise sub- 
stantially the temperature atmosphere before integrating the full 
equations. This is done during a time limited to f < 1: we inte- 
grate only the temperature equation (either with the full energy 
terms of the plain relaxation term, depending on the model), and 
the equation of hydrostatic equilibrium to obtain the pressure. 
We then use the equation of state to deduce the density from 
pressure and temperature. 

After this short hydrostatic heating phase, we integrate the 
full equations. As stated above, we also perform a depression 
at the outlet, which forms the transonic wind after some time. 
When using the full heat sources and sinks, we use successively 
two sets of parameters (see runs in Table 1). The aim is to ap- 
proach a density close to n = lO^cm"^ where the temperature 
reaches T = Q.5MK. The first set of parameters (run SI) leads 
to a density a bit below this value (see fourth panel in fig [3]i; 
the second set (all other runs) achieves the desired result. Note 
that the first choice of r;,„, - 10"* (run SI) limits somewhat the 
conductive flux, but that the second value (other runs) leads to no 
limitation at all. Figure[3]shows several profiles defining the final 
wind (with set number 2), as well as (bottom right panel) the two 
temperature/density curves obtained with both sets. A transition 
region is seen to occur at about 0.01 Rs - 7 000 km. The tem- 
perature has a minimum about 6000 K at the bottom, and a max- 
imum around 1.6 MK, around 4 solar radii. The chromospheric 
temperature curve is a smooth ramp, leading to a temperature 
close to 10000 K at the foot of the TR. This ramp is due to the 
artificial diff'usive term included in the temperature and pressure 
equations (3-4-5). The sonic Mach number is reached at about 5 
solar radii. Note that the heat sources and sinks balance (bottom 
left panel): conductive heating and radiative loss balance at the 
TR, while conductive loss and mechanical heating balance (with 
the help of the adiabatic cooling, not shown) in the corona. 



3.3. Wave transfer in the full wind model 

We now inject waves at the base of the previous stationary at- 
mosphere. We start from the relaxed wind (run S3) with density 




0.0001 0.0010 0.0100 0.1000 1.0000 10.0000 




Fig. 3. Stationary transonic wind with full heating sources and 
sinks (end of run S2, see Table 2). Top left: density and tem- 
perature; top right: velocity; bottom left: mechanical heating 
source (thick line), radiative cooling term (dotted), conductive 
term (dashed line); note sources and sinks are as they appear in 
the temperature equation, bottom right: density versus tempera- 
ture, run S2 with thick line, run SI with dotted line. 



about lO^cm ^ where the plasma has temperature T - 0.5 MK 
(see figure [3]). 

We consider as in Section 3.1 monochromatic waves: the ve- 
locity is given a sinusoidal pattern, after a transient of a half pe- 
riod. We deal only with the evanescent case, since it is the case 
which shows the largest diff'erence with the atmosphere without 
heat sources. The base period is fixed to be t = 10 min. Since 
the cut-off period below the TR is tq = 6 - 8 min, this implies 
that the wave is evanescent everywhere in the dense layers. The 
base amplitudes are either 4m/s or 40m/s. Top panels in Figure|4] 
show the low amplitude, the bottom panels the high amplitude. 

In both runs, the TR adopts as in the relaxation model a pe- 
riod of 6 min., about the cut-off period (fig.|5] The dotted fluc- 
tuations (not at scale) show the base velocity fluctuations in both 
runs, for comparison. 

Note that a peculiar behaviour of density and temperature is 
observed at the base: a transient systematic drift is superposed to 
the fluctuations. This drift modifies the mean gradients near the 
base, and allows the atmosphere to solve the contradiction im- 
posed by the incoming characteristics which have in-phase fluc- 
tuations of velocity, temperature, density, while upward propa- 
gating eigenmodes have phase-lags between the diff'erent quan- 
titites. 

Run B4 shows two phases, one with about three large os- 
cillations, while the rest of the fluctuations is of much smaller 
amplitude. Run B40 - with the largest amplitude - stops before 
the before second phase begins, because the profiles become too 
stiff for the numerical setup after the last temperature minimum, 
that is, after 60.25. (Although the calculation still proceeds up 
to time t=60.3, the time interval between 60.25 and 60.3 is un- 
physical, being invaded by numerical noise). Note the prominent 
nonlinear features of the large amplitude run: sawtooth velocity 
pattern, and peaks for the temperatures maxima, large tempera- 
ture troughs. 
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run B4 run 34 




0.001 0.010 0.100 0.00' O.QIO 0.100 

R - 1 R 1 

Fig. 4. Evanescent waves propagating through the full wind 
model. Top: run B4, base amplitude 4 m/s (time interval 
(60,60.5)); bottom: run B40, base amplitude 40 m/s (time in- 
terval (60,60.25); left: velocity profiles ponderated by density 
and sound speed ( y[{n)ucs)\ right: temperature profiles. Note run 
B40 stops after time t - 60.3 (see next figure). 



3.4. Relaxation model again 

To see whether the result is generic or not, we replace now the 
energy equation by the plain temperature equation with a cou- 
pling to a given temperature profile, as already considered in 
section 3.1: in other words, we return to the relaxation model. 
However, we start now from run S3, with the temperature pro- 
file obtained through the full energy equation. Hence, the differ- 
ence with the runs just discussed in the previous subsection is 
that only the adiabatic terms will react to the impinging wave, 
since the energy terms are not present. But, again, the starting 
TR structure will be the same. 

Recall that the extra coupling term in the temperature equa- 
tion is of the form -k{T - T") (eq. |7]). We insist that /c is a 
small parameter, which, in practice, gives this term a value which 
is never larger than 20% of the physical conductivity, hence it 
should not modify substantially the growth rate of the instabil- 
ity, if it is there. 

The result, as shown in top of Figure |6] shows that now the 
large amplitude wave (same as B40) arriving at the TR remain 
stable (run C40, top panels): the TR oscillations are actually 
much smaller than previously. Most importantly, when we in- 
crease the amplitude again by a factor 10 (run C400), we still 
obtain TR stability, with shocks propagating in the high corona 
(figure|6] bottom panels). 




Fig. 5. Evanescent waves propagating through the full wind 
model. Temperature and velocity at TR versus time; dotted: base 
velocity; top: run B4; bottom: run B40 



Fig. 6. Evanescent large-amplitude waves without energy terms, 
coupling with stationary temperature profile. Top: run C40, base 
amplitude 40 m/s; bottom: run C400, base amplitude 400 m/s 
left: velocity profiles ponderated by density and sound speed, 
y[in)ucs; right: temperature profiles. 



The amplitude being larger, one might suspect that the spa- 
tial scheme is not able to cope with too large and/or steep ve- 
locity/temperature/density profiles. However, a point deserves 
attention, namely the nonlinear relation between the oscillation 
amplitude at the TR and at the base: the strongest maximum of 
the velocity in run B4 is about 1.5 km/s, while the correspond- 
ing maximum is about 25 km/s, which is larger than the factor 
10 which relates the base amplitudes. Moreover, the beginning 
of the first phase is superficially similar to a (short) exponen- 
tial growth in both runs. So the hypothesis that the phenomenon 
we see is a true physical instability deserves to be examined in 
detail. 



4. Discussion 

It is rather clear the response of the TR to impinging transient 
evanescent waves is very different, depending on whether extra- 
adiabatic terms are present or not. The instability observed is 
nonlinear, that is, it shows up apparently only when waves have 
a sufficiently high amplitude. 

To establish more firmly the physical status of the instability, 
we redo the experiment with large amplitude wave (B40), but 
now with double resolution, dividing each mesh in two parts. 
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run D40 



run D40 



0.04 
0.03 
0.02 



0.00 
7( 




.0 70,1 70.2 70,3 70.4 70.5 
time 



run D80 



70,0 70,1 70,2 70,3 70,4 70,5 
time 



run D80 



0.04 

0.03 
0.02 
0.01 

0.00 
70 





70.1 70.2 70.3 70.4 70.5 
time 



70.0 70.1 70.2 70.3 70.4 70.5 
time 



Fig. 7. Evanescent waves, full energy equation, double resolu- 
tion; temperature and velocity at TR versus time; top: run D40, 
base amplitude 40m/s; bottom: run D80, base amplitude 80m/s 






Fig. 8. Evanescent waves, full energy equation, double resolu- 
tion; successive profiles of temperature (top), velocity (middle) 
and cooling term as in eq. (3) and (6) (bottom), in the vicinity 
of the TR during time interval 70.125 and 70.15; left: run D40, 
base amplitude 40m/s; right: run D80, base amplitude 80m/s 



Asking for the same domain size, we thus end up with - 599 
grid points instead of N = 300. 

The result is as shown in top of figure|7j stable (run D40). We 
thus increase the base amplitude by a factor two (run D80), and 
now recover the instability of velocity and temperature (bottom 
figure). It is interesting to examine in detail the two successive 



phases in run D40 (the stable one), because they show all char- 
acteristics of run B40 (nonlinearities), but during a longer time, 
since it escapes numerical catastrophy. During the first (nonlin- 
ear) phase, temperature and velocity are very different: tempera- 
ture shows symmetric spikes (only for maxima), while velocity 
shows sawtooths. During the second phase (which does'nt ex- 
ist for run D80), one finds ordinary quasi-monochromatic oscil- 
lations of smaller amplitude for both fields. The early velocity 
sawtooth pattern reveals shock formation, while the spiky ap- 
pearance of the temperature pattern reveals the (nonlinear) ad- 
vection of the temperature "wall" by the velocity field. There is 
a last (crucial) feature shown by temperature, that is, the dan- 
gerous decrease of the temperature troughs, which disappears in 
the second phase for run D40, and finally leads to the end of run 
D80. 

In fact, it seems that B40 is very close to run D40, the only 
difference being that D40 is able to pass the dangerous phase 
with the largest fluctuations, and not B40. So, we can safely con- 
clude that, while this difference is due to resolution change (and 
associated filtering change), the early nonlinear large response 
of the TR to upward propagating waves is physical, i.e., not due 
to the details of the spatial schema. 

Note that the observed instability has an amplitude threshold, 
but that the response probably depends on the incoming pertur- 
bation being transient, as the spectrum is changing with time as 
time elapses after starting the injection (the signal is not purely 
monochromatic, beginning at a finite time). Note also that the 
most striking difference between the low amplitude run (run B4) 
and the higher amplitude run (B40), as shown in figure |4] and 
also figure [5] is the appearance of large cold troughs at the foot 
oftheTR. 

While the temperature peaks (see figure |5] or figure |7]l can 
be attributed to the (nonlinear) advection of the TR temperature 
"wall" by the large velocity fluctuations, what is the origin of 
the (growing, unstable) temperature troughs? We propose that 
the cause both of the growth of the cold troughs is the radia- 
tive cooling term. This has the immediate advantage of explain- 
ing the nonlinearity of the response to the wave base amplitude. 
Indeed, the cooling term dependence on temperature is very non- 
linear in the low temperature range (see section 2), with a cut-off 
at IQ'^K (eq[l0|l. 

Here is the detailed scenario we propose. Consider the up- 
ward propagating wavefront resulting from a perturbation at the 
base, (for instance the beginning of our monochromatic injec- 
tion) which is steepening progressively at it propagates, at a rate 
increasing with altitude (as its amplitude is growing with height, 
due to rough energy density conservation). We know that this 
progressive wave is not really equivalent to a plain sound wave, 
since it is a gravito-acoustic one, but we will take this into ac- 
count later. 

If the wave was submitted only to the usual adiabatic pro- 
cesses, then it would proceed as just said: propagating, growing, 
steepening. But there are extra terms in the energy equation. The 
conductive term is definitely negligible below the TR, as well as 
the mechanical flux term. Note that the radiative cooling term 
remains strictly zero, but only as l ong as temperature remains 
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below the threshold (here IQ'^K, eq. 

Hence, this remains so all the way up to the TR, if the ini- 
tial amplitude of the wave is small enough for the threshold 
r* = lO'^K not to be trespassed during propagation. But, in the 
opposite case, then the wave crest is cooled immediately, with 
a rate proportional to the local density of the plasma. Hence the 
temperature maxima of a large amplitude wave will be cooled 
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both sooner (at lower altitude) and more rapidly than a small 
amplitude one. 



run D80 run D80 




0.01 0.01 
R-1 R-1 

Fig. 9. run D80: Instantaneous profiles of velocity, density, tem- 
perature and pressure at t=70.15; left: velocity (thick line), den- 
sity fluctuation (dotted), temperature fluctuation (dashed line); 
right: velocity (thick line), pressure fluctuation (dotted) 



What is the effect of the temperature drop when the wave 
crest is cooled down? Figure [8] shows for runs D40 et D80 
10 successive profiles at early times, from time t=70.125 to 
t=70. 15, of temperature, velocity and cooling term in the vicinity 
of the TR. The time interval corresponds to the first large oscil- 
lation at the TR, as seen in figure |7] It is seen that the growth of 
the cooling term and the growth of the temperature and velocity 
wavefront go together. Recall that both the conductive and me- 
chanical heating terms are negligible in the region discussed, as 
stated above. 

Now, consider for a moment that gravity can be neglected, 
i.e., that the wave frequency is high. Assume again local tem- 
perature is just below the threshold for radiative cooling. Also, 
consider as initial conditions not a progressive wave, but, in- 
stead, a local compressive movement (with no initial pressure 
fluctuation). This movement indeed will be accelerated as soon 
as the cooling threshold is trespassed, because the pressure in- 
crease which should resist the compressive movement is sud- 
denly missing. 

But these conditions are not those of the progressive waves 
propagating from the base. This is because (If gravity is ne- 
glected) then pressure, temperature, density and velocity fluc- 
tuations are in phase; in that case, cooling of the temperature 
maxima will affect the wave, but will not lead to an unstable 
compression. 

What finally explains the instability in the present problem 
is that we are dealing not with a high-frequency wave, but with a 
wave at frequency close to the cut-off. In that case, density, pres- 
sure, velocity are no longer in phase, as may be seen in figure |9] 
which shows in details the different profiles at the last time of the 
series shown in figure [8] The left panel shows that the temper- 
ature crest is leading the velocity crest, while the density crest 
lies behind it. If we imagine now that the temperature crest is 
decreased by the radiative cooling, this will decrease the right 
foot of the pressure bump which should maintain the wave form. 
Hence the fluid ahead of the wave is accelerated toward the pres- 
sure drop, which is exactly what we see. This in turn raises the 
pressure, hence pushes again the temperature above the thresh- 
old, etc... 

This argument explains also why injecting evanescent waves 
leads to instability, and not injecting propagative waves. This is 
because injecting evanescent waves generate much more waves 
with the cut-off frequency than the propagative waves (Fleck and 



Schmitz 1991, Kalkofen et al 1994) and because the instability 
occurring here requires that temperature, density and velocity be 
out of phase as found above in the transient propagative waves 
at the cut-off frequency. 

Now, this instability clearly stops working when it enters the 
TR. This is because In the TR, conductive terms become impor- 
tant, and are able to balance cooling, so the instability stops. 

Hence, the instability has nothing to do with the TR it- 
self: it is not an instability of the TR, but the instability of 
the (quasi-)isothermal, low temperature part of the atmosphere, 
which should occur whenever there is a layer which is optically 
thin, and cold enough, so that radiative cooling is active, and 
only it. Finally, the TR stops the instability, because of the stabi- 
lizing effect of thermal conduction there. 

Also, in the cases (say, D40) where the evolution can be fol- 
lowed all the through the unstable phase, why does the process 
not repeat indefinitely for the next wavefronts arriving from the 
base? Why does the unstable phase stop after two or three peri- 
ods? We think this last point has to do with the progressive dis- 
appearance of the cut-off part of the spectrum with time (Fleck 
and Schmitz (1991), Kalkofen et al. (1994)), that is, it is not, 
again, related to the TR itself. 

The present scenario explains well our simulations, but what 
about the real sun? In the real sun, there is a physical limit to the 
radiative cooling term which is used here, which has not only 
a temperature threshold as considered here, but also a density 
threshold, since in the densest layers the plasma becomes opti- 
cally thick. 

There are several ways to take this into account. The simplest 
is to introduce a dependence of the cooling term on density in the 
pressure equation. A first possibility is as follows: -n^MJ) be- 
coming n^n* l{n + n*)A{T). With this prescription, the cooling 
rate becomes density-independent when density is larger than 
n* . Note that if we set n* close to the TR density itself (typi- 
cally, 10** COT"''), its very structure will be strongly modified. We 
have chosen n* = 10"cot"'', and found that the growth rate of 
the first oscillations in a run equivalent do D80 was not basi- 
cally changed. A further case would be to replace -n^A{T) by 
(nn* /{n + n*)^A(T), which would completely switch off the term 
at layers denser than n*. 

To summarize, the upward propagating wave packets gen- 
erated by base perturbations will implode only in the layers 
which are comprized between the optically thick region and the 
TR where conduction becomes large, which stops the instabil- 
ity. Hence the instability should be most observable when and 
where the TR is high enough above the optically thick layers. 
This should provide an observational check of the present work. 

Work is in progress to include magnetic effects, and general- 
ize our study to an axisymmetric corona. 
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